library(tidyverse)
library(data.table)
library(readxl)
library(GGally)
library(corrplot)
library(leaflet)
library(raster)
library(mltools)
library(highcharter)
raw_encuesta <- fread("databases/calidad_vida_ok.csv", encoding = "UTF-8") %>%
as_tibble()
# Unificación de las estadísticas de los barrios "ALTAVISTA" y "CABECERA ALTAVISTA"
raw_encuesta <- raw_encuesta %>% mutate(
encuesta_calidad.barrio = case_when(
encuesta_calidad.barrio == "CABECERA ALTAVISTA" ~ "ALTAVISTA",
encuesta_calidad.barrio == "CIUDADELA NUEVO OCCIDENTE" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
encuesta_calidad.barrio == "AREA DE EXPANCION SAN CRISTOBAL" ~ "ÁREA DE EXPANSIÓN SAN CRISTÓBAL",
encuesta_calidad.barrio == "CABECERA SAN CRISTÓBAL" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
encuesta_calidad.barrio == "SAN CRISTOBAL" ~ "CABECERA URBANA CORREGIMIENTO SAN CRISTÓBAL",
encuesta_calidad.barrio == "PROGRESO Nº 2" ~ "PROGRESO",
TRUE ~ encuesta_calidad.barrio
)
)
diccionario <- read_excel("databases/diccionario_ecv_2016.xlsx")
raw_encuesta %>%
count(`encuesta_calidad.barrio`, sort = T) %>%
arrange(n)
## # A tibble: 296 x 2
## encuesta_calidad.barrio n
## <chr> <int>
## 1 EL CARMELO 4
## 2 GUAYAQUIL 4
## 3 LA ILUSIÓN 6
## 4 SAN JOSÉ 10
## 5 BOQUERÓN 15
## 6 EL PATIO 19
## 7 BARRIO COLOMBIA 21
## 8 DESCONOCIDO 21
## 9 YARUMALITO 22
## 10 EL UVITO 25
## # ... with 286 more rows
encuesta_ingresos<-raw_encuesta%>%
mutate(proporcion_ingreso_extra=case_when(encuesta_calidad.p_71 == 1 ~1,
TRUE~ 0))%>%
mutate(proporcion_pension=case_when(encuesta_calidad.p_71 != 0~1,
TRUE~ 0))%>%
mutate(subsidio=case_when(encuesta_calidad.p_93 != 0~1,
TRUE~ 0))%>%
mutate(salario=case_when(encuesta_calidad.p_121 != 0~1,
TRUE~ 0))%>%
mutate(arriendo=case_when(encuesta_calidad.p_227 != 0~1,
TRUE~ 0))%>%
mutate(total_gastos=case_when(encuesta_calidad.p_232 != 0~1,
TRUE~ 0))%>%
mutate(gasto_servicios=case_when(encuesta_calidad.p_240 != 0~1,
TRUE~ 0))%>%
mutate(gasto_hogar=case_when(encuesta_calidad.p_254 != 0~1,
TRUE~ 0))%>%
mutate(encuesta_calidad.barrio = str_replace(encuesta_calidad.barrio, "ANDALUCIA", "ANDALUCÍA") %>%
str_replace("Nº 2", "NO.2") %>%
str_replace("Nº 1", "NO.1") %>%
str_replace("Nº 3", "NO.3") %>%
str_replace("AREA EXPANSION", "ÁREA DE EXPANSIÓN") %>%
str_replace("EXPANCION", "EXPANSIÓN") %>%
str_replace("AREA", "ÁREA") %>%
str_replace("BOMBONA", "BOMBONÁ") %>%
str_replace("LA ASOMADERA", "ASOMADERA") %>%
str_replace("BELALCAZAR", "BELALCÁZAR") %>%
str_replace("CALAZANS", "CALASANZ") %>%
str_replace("COLON", "COLÓN") %>%
str_replace("MIRA FLORES", "MIRAFLORES") %>%
str_replace("BARRIO FACULTAD DE MINAS", "FACULTAD DE MINAS") %>%
str_replace("CABECERA SAN ANT DE PR.", "SAN ANTONIO DE PRADO") %>%
str_replace("CARLOS E RESTREPO", "CARLOS E. RESTREPO") %>%
str_replace("URQUITA", "URQUITÁ") %>%
str_replace("LOS CERROS EL VERJEL", "LOS CERROS EL VERGEL") %>%
str_replace("CAYCEDO", "CAICEDO") %>%
str_replace("VALDES", "VALDÉS") %>%
str_replace("CERRO EL VOLADOR", "B. CERRO EL VOLADOR") %>%
str_replace("MOSCU", "MOSCÚ") %>%
str_replace("JOSELA", "JOSÉ LA") %>%
str_replace("JOSE", "JOSÉ") %>%
str_replace("EL YOLOMBO", "YOLOMBO") %>%
str_replace("PIEDRAS BLANCAS", "PIEDRAS BLANCAS - MATASANO") %>%
str_replace("BASILIA", "BRASILIA") %>%
str_replace("VILLA TINA", "VILLATINA") %>%
str_replace("LILIAM", "LILLIAM") %>%
str_replace("BOLIVAR", "BOLÍVAR") %>%
str_replace("CORREGIMIENTO PALMITAS", "PALMITAS SECTOR CENTRAL") %>%
str_replace("INES", "INÉS") %>%
str_replace("FE", "FÉ") %>%
str_replace("LUCIA", "LUCÍA") %>%
str_replace("SABIO", "SAVIO") %>%
str_replace("BERMEJAL- LOS ÁLAMOS", "BERMEJAL-LOS ÁLAMOS") %>%
str_replace("BOLÍVARIANA", "BOLIVARIANA") %>%
str_replace("EL NOGAL - LOS ALMENDROS", "EL NOGAL-LOS ALMENDROS") %>%
str_replace("JUAN XXIII - LA QUIEBRA", "JUAN XXIII LA QUIEBRA") %>%
str_replace("PROGRESO Nº 2", "EL PROGRESO") %>%
str_replace("MARIA", "MARÍA") %>%
str_replace("PLAYÓN", "PLAYON") %>%
str_replace("EL SOCORRO / LA GABRIELA", "EL SOCORRO") %>%
str_replace("FÉRRINI", "FERRINI") %>%
str_replace("LA CANDE LARIA", "LA CANDELARIA") %>%
str_replace("EL PLAYON", "PLAYÓN") %>%
str_replace("IGUANA", "IGUANÁ") %>%
str_replace("MARÍA CANO - CARAMBOLAS", "MARÍA CANO-CARAMBOLAS") %>%
str_replace("DE ABURRA", "DEL ABURRÁ") %>%
str_replace("ALTAVISTA CENTRAL", "ALTAVISTA SECTOR CENTRAL") %>%
str_replace("SECTOR CENTRAL", "CENTRO ADMINISTRATIVO") %>%
str_replace("ALTAVISTA CENTRO ADMINISTRATIVO", "ALTAVISTA SECTOR CENTRAL") %>%
str_replace("SANTA ELENA CENTRO ADMINISTRATIVO", "SANTA ELENA SECTOR CENTRAL") %>%
str_replace("PALMITAS CENTRO ADMINISTRATIVO", "PALMITAS SECTOR CENTRAL") %>%
str_replace("PROGRESO", "EL PROGRESO")
)
db_ingresos<-encuesta_ingresos%>%
group_by(encuesta_calidad.barrio, encuesta_calidad.comuna)%>%
summarize(n=n(),
proporcion_ingreso_extra=sum(proporcion_ingreso_extra==1,na.rm=TRUE)/sum(!is.na(proporcion_ingreso_extra),na.rm=TRUE),
proporcion_pension=sum(proporcion_pension==1,na.rm = TRUE)/sum(!is.na(proporcion_pension),na.rm=TRUE),
subsidio=sum(`encuesta_calidad.p_93`, na.rm=TRUE)/sum(subsidio==1,na.rm = TRUE),
salario=sum(`encuesta_calidad.p_121`, na.rm=TRUE)/sum(salario==1,na.rm = TRUE),
arriendo=sum(`encuesta_calidad.p_227`,na.rm=TRUE)/sum(arriendo==1,na.rm = TRUE),
total_gastos=sum(`encuesta_calidad.p_232`,na.rm=TRUE)/sum(total_gastos==1,na.rm = TRUE),
gasto_servicios=sum(`encuesta_calidad.p_240`,na.rm=TRUE)/sum(gasto_servicios==1,na.rm = TRUE)
)%>%ungroup()%>%
mutate(subsidio= replace_na(subsidio, 0))%>%
mutate(salario= replace_na(salario, 0))%>%
mutate(arriendo= replace_na(arriendo, 0))%>%
mutate(total_gastos= replace_na(total_gastos, 0))%>%
mutate(gasto_servicios= replace_na(gasto_servicios, 0))
# Eliminamos el barrio Desconocido
db_ingresos <- db_ingresos %>%
filter(encuesta_calidad.barrio != "DESCONOCIDO")
db_ingresos[complete.cases(db_ingresos),]
## # A tibble: 308 x 10
## encuesta_calida~ encuesta_calida~ n proporcion_ingr~
## <chr> <chr> <int> <dbl>
## 1 AGUAS FRÍAS ALTAVISTA 126 0.00794
## 2 ALDEA PABLO VI POPULAR 426 0.00235
## 3 ALEJANDRÍA EL POBLADO 709 0.00987
## 4 ALEJANDRO ECHAV~ BUENOS AIRES 1094 0.00640
## 5 ALFONSO LÓPEZ CASTILLA 1994 0.00752
## 6 ALTAMIRA ROBLEDO 755 0.00927
## 7 ALTAVISTA ALTAVISTA 168 0.0119
## 8 ALTAVISTA BELEN 1988 0.00553
## 9 ALTAVISTA SECTO~ ALTAVISTA 507 0.00592
## 10 ALTOS DEL POBLA~ EL POBLADO 375 0.00267
## # ... with 298 more rows, and 6 more variables: proporcion_pension <dbl>,
## # subsidio <dbl>, salario <dbl>, arriendo <dbl>, total_gastos <dbl>,
## # gasto_servicios <dbl>
Escritura.
write_excel_csv2(db_ingresos, "databases/db_ingresos.csv")
write_excel_csv2(db_ingresos[complete.cases(db_ingresos),], "databases/db_ingresos_completecases.csv")
Esta dimensión llamada “ingresos” aborda algunas preguntas referentes a los ingresos y gastos de las personas de los diferentes barrios de Medellín.
gasto_servicio:Es el promedio por cada barrio de la cantidad del ingreso mensual destina a servicios públicos en cada hogar.
Nota:En la base de datos utilizada se encontraron algunas celdas vacías en las preguntas utilizadas para calcular las variables “subsidio”, “salario”, y “arriendo”, esto puede deberse a que la persona encuestada no recibe algún tipo de subsidio o salario, y no paga arriendo. Al fijarnos en los barrios en los que para ninguna persona hay respuesta, nos dimos cuenta que son barrios de invasión, por lo cual es coherente que haya un cero para tal respuesta; también se pudo observar que en algunos casos las personas no reciben sueldo, pero sí reciben un subsidio. Por ésto, en las celdas de nuestra nueva base de datos db_ingresos en donde no fue posible calcular promedios, pues el cociente quedaba en cero y esto tiende a infinito, se reemplazó por ceros.
summary(db_ingresos)
## encuesta_calidad.barrio encuesta_calidad.comuna n
## Length:308 Length:308 Min. : 4.0
## Class :character Class :character 1st Qu.: 348.8
## Mode :character Mode :character Median : 862.5
## Mean :1073.2
## 3rd Qu.:1439.0
## Max. :8987.0
## proporcion_ingreso_extra proporcion_pension subsidio
## Min. :0.000000 Min. :0.1000 Min. : 0
## 1st Qu.:0.003061 1st Qu.:0.4049 1st Qu.: 244619
## Median :0.005511 Median :0.4426 Median : 370100
## Mean :0.006463 Mean :0.4386 Mean : 520659
## 3rd Qu.:0.008005 3rd Qu.:0.4696 3rd Qu.: 587997
## Max. :0.066667 Max. :0.7500 Max. :8241967
## salario arriendo total_gastos gasto_servicios
## Min. : 0 Min. : 0 Min. : 425763 Min. : 50500
## 1st Qu.: 305632 1st Qu.: 283998 1st Qu.: 833227 1st Qu.:150364
## Median : 437404 Median : 409747 Median :1087169 Median :206197
## Mean : 703446 Mean : 535185 Mean :1393247 Mean :225134
## 3rd Qu.: 742865 3rd Qu.: 654791 3rd Qu.:1580345 3rd Qu.:275162
## Max. :8734250 Max. :1901874 Max. :4976319 Max. :588615
min_gasto_servicios<-apply(db_ingresos[10],2,min, na.rm=T)
db_ingresos[which(db_ingresos$gasto_servicios==min_gasto_servicios),1]
## # A tibble: 1 x 1
## encuesta_calidad.barrio
## <chr>
## 1 SAN JOSÉ
max_gasto_servicios<-apply(db_ingresos[10],2,max, na.rm=T)
db_ingresos[which(db_ingresos$gasto_servicios==max_gasto_servicios),1]
## # A tibble: 1 x 1
## encuesta_calidad.barrio
## <chr>
## 1 ALTOS DEL POBLADO
grafico<-ggpairs(db_ingresos,columns = 4:10, mapping=aes(alpha = 0.4))
ggsave(filename="images/correlacion.png", plot =grafico, width = 20, height = 20)
Correlación entre las variables
data<-as.matrix(db_ingresos[,3:9])
mcor<-cor(data)
mcor
## n proporcion_ingreso_extra
## n 1.00000000 -0.06518933
## proporcion_ingreso_extra -0.06518933 1.00000000
## proporcion_pension -0.17011586 -0.03081172
## subsidio -0.01152908 -0.02913350
## salario -0.02949656 -0.01847511
## arriendo -0.10414827 0.01745026
## total_gastos -0.10706213 0.01250525
## proporcion_pension subsidio salario
## n -0.17011586 -0.01152908 -0.02949656
## proporcion_ingreso_extra -0.03081172 -0.02913350 -0.01847511
## proporcion_pension 1.00000000 0.03797035 0.08478863
## subsidio 0.03797035 1.00000000 0.28460736
## salario 0.08478863 0.28460736 1.00000000
## arriendo 0.05240449 0.40708611 0.65315248
## total_gastos 0.02899374 0.40252794 0.63576569
## arriendo total_gastos
## n -0.10414827 -0.10706213
## proporcion_ingreso_extra 0.01745026 0.01250525
## proporcion_pension 0.05240449 0.02899374
## subsidio 0.40708611 0.40252794
## salario 0.65315248 0.63576569
## arriendo 1.00000000 0.95678984
## total_gastos 0.95678984 1.00000000
corrplot(mcor,type = "upper",order = "hclust",tl.col = "black",tl.srt = 45)
db_ingresos_escalado <- db_ingresos%>%
mutate_at(vars(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna), scale)
datos_escalados_agrupamiento <- db_ingresos_escalado %>%
dplyr::select(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna)
pca <- prcomp(datos_escalados_agrupamiento)
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.8902 1.0168 0.9845 0.8853 0.68540 0.3579 0.20592
## Proportion of Variance 0.5104 0.1477 0.1385 0.1120 0.06711 0.0183 0.00606
## Cumulative Proportion 0.5104 0.6581 0.7966 0.9085 0.97564 0.9939 1.00000
porcentaje_exp <- 100 * (pca$sdev ^ 2) / sum(pca$sdev ^ 2)
ggplot(data = data.frame(componente = factor(1:ncol(datos_escalados_agrupamiento)), porcentaje_exp)) +
geom_col(mapping = aes(x = componente, y = porcentaje_exp))
suma_acumulada_porcentaje <- cumsum(porcentaje_exp)
ggplot(data = data.frame(componente = factor(1:ncol(datos_escalados_agrupamiento)), suma_acumulada_porcentaje),
aes(x = componente, y = suma_acumulada_porcentaje, group = 1)) +
geom_point() +
geom_line() +
geom_hline(yintercept = 90, color = "red") +
geom_vline(xintercept = 7)
suma_acumulada_porcentaje
## [1] 51.03878 65.80981 79.65645 90.85302 97.56407 99.39425 100.00000
datos_proyectados <- t(t(pca$rotation) %*% t(datos_escalados_agrupamiento))
datos_agrupamiento <- datos_proyectados[, 1:4] %>% as_tibble()
Se utilizará la metodología de Elbow para identificar el k óptimo que ofreza una menor distancia entre clusters como medida de error Wss.
datos_escalados_agrupamiento<- db_ingresos_escalado %>%
dplyr::select(-n, -encuesta_calidad.barrio, -encuesta_calidad.comuna)
k_maximo <- 15
wss <- sapply(1:k_maximo,
function(k) {
kmeans(datos_escalados_agrupamiento, k, nstart = 10, iter.max = 15 )$tot.withinss
})
plot(1:k_maximo,wss,
type="b", pch = 19, frame = FALSE,
xlab="Número de grupos [K]",
ylab="Total de suma de cuadrados intra-grupos")
Según nuestro algoritmo de K-means, nuestro K óptimo es de 9 grupos, pero se considera más pertinente ultilizar un número menor para favorecer la interpretabilidad.
k_optimo_kmeans<-9
k_seleccionado<-5
Fijaremos una semilla de los número aleatorias para garantizar la reproducibilidad de las conclusiones del agrupamiento.
set.seed(3)
agrupamiento_kmeans <- kmeans(datos_agrupamiento, centers = k_seleccionado, nstart = 10, iter.max = 15 )
db_ingresos_kmeans <- bind_cols(db_ingresos, cluster = agrupamiento_kmeans$cluster)
political <- shapefile("Barrio_Vereda/Barrio_Vereda.shp")
Encoding(political@data$NOMBRE) <- "UTF-8"
political$NOMBRE <- political$NOMBRE %>% toupper() %>% str_replace("DE MESA", "DE MESA")
grupos_barrios <- data.frame(barrio_nombre = db_ingresos$encuesta_calidad.barrio, grupo = agrupamiento_kmeans$cluster)
nombres_mapa <- data.frame(nombre_barrio = political$NOMBRE)
vector_nombres = c()
vector_grupos = c()
for(nombre_mapa in nombres_mapa$nombre_barrio) {
grupo <- grupos_barrios[grupos_barrios$barrio_nombre == nombre_mapa, 2][1]
vector_nombres <- c(vector_nombres, nombre_mapa)
vector_grupos <- c(vector_grupos, grupo)
}
factpal <- colorFactor(rainbow(k_seleccionado), vector_grupos)
mapa_grupos <-tibble(
nombre_barrio = vector_nombres,
grupo = vector_grupos,
color = factpal(grupo),
descripcion = case_when(
grupo == 1 ~ "Grupo 1",
grupo == 2 ~ "Grupo 2",
grupo == 3 ~ "Grupo 3",
grupo == 4 ~ "Grupo 4",
grupo == 5 ~ "Grupo 5",
TRUE ~ "Sin muestra representativa")
)
colores <- mapa_grupos %>%
dplyr::select(-nombre_barrio) %>%
distinct(color, .keep_all = TRUE) %>%
arrange(desc(grupo))
leaflet(data = political) %>%
addTiles() %>%
addProviderTiles(providers$OpenStreetMap) %>%
addPolygons(fill = TRUE, stroke = TRUE, weight = 2, color = mapa_grupos$color,
label = as.character(political$NOMBRE),
popup = as.character(political$NOMBRE)) %>%
addLegend("bottomright", colors = colores$color, labels = as.character(colores$descripcion))